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ABSTRACT 

The modeling of galaxy formation and reionization, two central issues of modern cosmology, relies 
on the accurate follow-up of the intergalactic medium (IGM). Unfortunately, owing to the complex 
nature of this medium, the differential equations governing its ionization state and temperature are 
only approximate. In this paper, we improve these master equations. We derive new expressions for 
the distinct composite inhomogeneous IGM phases, including all relevant ionizing/recombining and 
cooling/heating mechanisms, taking into account inflows/outflows into/from halos, and using more 
accurate recombination coefficients. Furthermore, to better compute the source functions in the equa¬ 
tions we provide an analytic procedure for calculating the halo mass function in ionized environments, 
accounting for the bias due to the ionization state of their environment. Such an improved treatment 
of IGM evolution is part of a complete realistic model of galaxy formation presented elsewhere. 
Subject headings: cosmology: theory — intergalactic medium — galaxies — galaxies: formation 


1. INTRODUCTION 

The evolution of galaxies is intertwined with that of 
the intergalactic medium (IGM). Mechanical heating of 
IGM by active galactic nuclei (Bower et al. 2006; Croton 
et al. 2006) and radiative heating by X-rays produced 
in supernovae (White & Rees 1978; Dekel & Silk 1986; 
Cole 1991; White & Frenk 1991; Lacey & Silk 1991; Oh 
& Haiman 2003) together with ionizing photons emitted 
by young stars (Ikeuchi 1986; Rees 1986; Shapiro et al. 
1990; Miralda-Escude & Ostriker 1990, 1992; Efstathiou 
1992) modify the temperature and ionization state of the 
IGM, which in turn alters subsequent galaxy formation. 

The physics involved in the coupled evolution of IGM 
and luminous sources is so complex and covers such a 
wide range of scales that its treatment involves important 
approximations. In fact, most studies focusing on galaxy 
formation adopt an IGM with fixed adhoc properties. 
Only studies of reionization do follow the IGM evolution 
in more or less detail. 

IGM evolution is described by a couple of differential 
equations for its ionization state and temperature with 
some source functions provided by a galaxy model. It is 
in this latter part where most approximations and simpli¬ 
fying assumptions are made, depending on the particu¬ 
lar approach followed, namely hydrodynamic simulations 
(Quinn, Katz, & Efstathiou 1996; Weinberg, Hernquist, 
& Katz 1997; Navarro & Steinmetz 1997; Ciardi et al. 
2000; Wyithe & Loeb 2003; Iliev et al. 2007; Okamoto et 
al. 2008; Trac et al. 2008; Battaglia et al. 2013; Sobacchi 
& Mesinger 2013a), numerical and seminumerical simu¬ 
lations (Zhang et al. 2007; Faucher-Giguere et al. 2009; 
Zahn et al. 2011; Sobacchi & Mesinger 2013b), pure an¬ 
alytic models (Haiman et al. 1996; Thoul & Weinberg 
1996; Dijkstra et al. 2004; Furlanetto et al. 2004; Alvarez 
et al. 2012; Kaurov & Gnedin 2013), and semianalytic 
models (Babul & Rees 1992; Efstathiou 1992; Shapiro 
et al. 1994; Mesinger & Dijkstra 2008; Font et al. 2011; 
Wyithe & Loeb 2013), each with its pros and cons. 

The treatment of the IGM itself, a composite inhomo¬ 
geneous multiphase medium, is not fully accurate either. 


In principle, the problem is less severe for hydrodynamic 
simulations than for (semi)numerical and (semi)analytic 
models because these equations apply locally, so one 
must not worry about the spatially fluctuating proper¬ 
ties of IGM. However, current simulations do not resolve 
the different ionized phases (Finlator et al. 2012). 

A usual procedure (e.g. Shapiro et al. 1994; Wyithe 
& Loeb 2003; Benson et al. 2006; Zhang et al. 2007) 
is to consider the IGM as having a simple hydrogenic 
composition and constant, uniform temperature, equal 
to the characteristic temperature of photoionized hydro¬ 
genic gas (~ 10 4 K), and to focus on the evolution of the 
ionization state through the simple equation derived by 
Shapiro & Giroux (1987). But the IGM temperature is 
crucial not only for estimating the minimum galaxy mass 
but also for computing the recombination coefficients, so 
such an approximation also affects the ionization state of 
the IGM. 

Hui & Gnedin (1997) derived the first coupled equa¬ 
tions for the ionization state and temperature of the IGM 
taking into account the dependence of the latter on hy¬ 
drogen and helium abundances and local density of the 
gas (Miralda-Escude & Rees 1994). However, these equa¬ 
tions only held for the cooling phase after ionization, and 
Haiman & Holder (2003) and Hui & Haiman (2003) ex¬ 
tended them to include the ionization period. 

But the IGM is also nmltiphasic (Miralda-Escude, 
Haehnelt, & Rees 2000 and references therein): the neu¬ 
tral, singly, and doubly ionized regions are separated. 
Choudhury & Ferrara (2005) derived the equations for 
IGM evolution since the dark ages taking into account 
the full composite, inhomogeneous, and multiphase na¬ 
ture of IGM. However, instead of taking the average re¬ 
combination coefficients in each (ionized) phase they use 
the value these coefficients would take for the average 
(approximately mass-weighted) IGM temperature. On 
the other hand, they ignored the mass exchanges be¬ 
tween halos and IGM, although about 90% of the initial 
diffuse gas ends up locked into halos, and the current 
IGM metallicity shows that halos also eject substantial 
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amounts of gas into the medium. 

These mass exchanges affect the volume filling factors 
of the various ionized species as well as the mean parti¬ 
cle kinetic energy, so they must be taken into account. 
In principle, this would introduce one explicit differen¬ 
tial equation for each of the varying comoving densities. 
But, taking into account their trivial form (i.e. the varia¬ 
tion in each quantity is equal to the corresponding source 
function), these variations can be directly included in the 
usual master equations. Note that the IGM metallic- 
ity determining its mean molecular weight also changes. 
However, the mass fraction in metals in the IGM is so 
small (~ 10 -2 Z 0 at z ~ 5; e.g. Sirncoe et al. 2011 
and D’Odorico et al. 2013) that these variations have a 
negligible effect. 

Lastly, the source functions in the IGM master equa¬ 
tions were calculated by averaging the feedback of lu¬ 
minous objects over ionized regions, assuming an evolv¬ 
ing universal halo mass function (MF). Yet, as the mass 
of halos able to trap gas and to form stars depends on 
the temperature and ionization state of the surround¬ 
ing IGM, the halo MF itself depends on the environ¬ 
ment. That is, the MF of halos lying in ionized or 
neutral regions differs. This bias, hereafter referred to 
as the ionization-bias to distinguish it from the well- 
known mass-bias (e.g. Tinker et al. 2010 and references 
therein) 1 must thus be corrected for. 

The aim of the present paper is to improve the analyt¬ 
ical treatment of IGM evolution by deriving new more 
accurate master equations for its ionization state and 
temperature, and by estimating the halo ionization-bias 
necessary to properly compute the source functions in 
these equations. Such an improved treatment of IGM 
can be incorporated into any given (semi(numerical or 
(semi(analytic model of galaxy formation such as the one 
developed by Manrique et al. (2015). The IGM prop¬ 
erties shown throughout the paper to illustrate the ef¬ 
fects of the new treatment have been obtained from that 
model. 

In Sections 2 and 3, we derive the new equations for the 
IGM ionization state and temperature, respectively. In 
Section 4, we derive the halo mass functions that result 
in neutral and ionized environments. Our results are 
discussed and summarized in Section 5. 

2. IONIZATION STATE EQUATIONS 

The structure of IGM is determined by the ionizing ra¬ 
diation from luminous sources. UV photons with a short 
mean free path ionize small regions around these sources. 
Their less energetic fraction gives rise to singly ionized 
hydrogen and helium bubbles, while the less abundant, 
more energetic fraction gives rise to doubly ionized he¬ 
lium subbubbles. Bubbles and subbubbles grow and pro¬ 
gressively overlap or retract and fragment, depending on 
the intensity of the ionizing flux is. In any case, the 
neutral, singly and doubly ionized phases are kept well 
separated at any time. 

As mentioned, IGM is not only multiphasic but also 
inhomogeneous. All IGM properties, such as tempera¬ 
ture, baryon density or Hi number density, are random 
fields characterized by their respective probability distri- 

1 The mass-bias is the dependence on large-scale mean density 
of the abundance of halos with a given mass. 


bution functions (PDFs). We are here interested in the 
time evolution of the IGM properties averaged over dif¬ 
ferent regions. When these averages refer to the neutral, 
singly, and doubly ionized phases, they will be denoted 
by angular brackets with subscripts I, II, and III, re¬ 
spectively; when they refer to regions encompassing one 
particular chemical species, such as Hii (i.e. all ionized 
regions), the subscript will explicitly indicate that chemi¬ 
cal species; and when the average is over the entire IGM, 
there will be no subscript. Averages of the product of 
several (either correlated or uncorrelated) quantities are 
for their joint PDF, so they will differ in general from 
the product of the averages of the individual quantities. 

The local comoving density of Hu ions, rimi, at the 
cosmic time t satisfies the balance equation 


dnmi 
d t 


Nmi 


am (T) 

-o— uhii n e , 

a A 


( 1 ) 


where n e is the comoving density of free electrons, IVhii is 
the local metagalactic emissivity of H i-ionizing photons 
due to luminous sources and recombinations, including 
redshifted photons emitted and not absorbed at higher 
z’s, and the second term on the right is the recombi¬ 
nation rate density to Hi. Note that the temperature- 
dependent recombination coefficient for optically thin re¬ 
gions, aHi(T’) (see e.g. Meiksin 2009; Fauclier-Giguere et 
al. 2009), is divided by the cube of the cosmic scale factor 
a(t) so as to express it in comoving units. 

Taking the average of equation (1) over the whole IGM, 
with the average of the second term on the right decom¬ 
posed in the sum of the averages over the different phases 
I, II, and III, duly weighted by their respective volume 
filling factors, Qi = 1 — Qhii, Qii = Qhii — QHein, and 
Qm = QHein, with Qhii and Qneiii standing for the 
Hu and Hem volume filling factors, respectively defined 
as (?rHii}/(nH) and (n H ein)/(nHe), gives rise to the rig¬ 
orous equation 
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Approximating o:hi(T) in ionized regions by a uniform 
value corresponding to the characteristic temperature 
Ttyp of photoionized hydrogenic gas (~ 10 4 K), and divid¬ 
ing by the approximately constant value (ignoring inflows 
and outflows) of the mean comoving hydrogen density, 
(uh), we arrive at the following simple equation for the 
Hu volume filling factor Qhii (Shapiro & Giroux 1987), 


dQffli (Nmi) a HI (T typ ) 

-r:— — ~I -\-5-U HII We/HII Q/ffll > Wl 

dr Wh) er 

where Chii = ( n mi)Hii/( n Hii)Hii 1S th e so-called clump¬ 
ing factor. To write equation (3), we have made two ap¬ 
proximations: (n e ) HII « (n H ii>Hii and (n H ) H n ~ («h)- 
The former presumes hydrogenic composition, and the 
latter presumes that ionized regions have the same aver¬ 
age properties as the whole IGM. 

However, (nn) is not constant, but evolves due to in¬ 
flows and outflows into and from halos. In addition, the 
IGM is not strictly hydrogenic, as its temperature varies 
both in space and time. Lastly, there should be, as men¬ 
tioned earlier, some halo ionization-bias, so the average 
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IGM properties in ionized regions should differ in general 
from the global average properties. We should thus try 
to do better. 

Let us comeback to the rigorous equation (2). Neglect¬ 
ing metals, we have n e = niiii(l + Y/4X) + riHeiii, where 
the comoving density of Hem ions, riHeim takes the ap¬ 
proximate form f(X, Y, 7 )nmi, with / equal to a univer¬ 
sal function of the hydrogen and helium mass fractions, 
X and Y, respectively, and the typical spectral index 7 
of ionizing sources. Thus, the average in the summation 
on the right of equation ( 2 ) splits into a sum of two prod¬ 
ucts of the form: average of a function of T times average 
of jj. This is possible thanks to the fact that there is 
essentially no correlation between T and tihii- The rea¬ 
son for this is that, in ionized regions, nnii is essentially 
equal to riu =In, where n is the baryon density. Fur¬ 
thermore, the only terms in equation (7) for the evolu¬ 
tion of the IGM temperature coupling n and T are the 
second and fifth ones giving the heating/cooling by adi¬ 
abatic compression/expansion of the fluid element, and 
the heating/cooling by the loss/gain of baryons due to in¬ 
flows/outflows, respectively, which are less than the first 
term giving the cosmic adiabatic cooling, and much less 
than the third and fourth terms including the stochastic 
effects of nearby luminous sources. Under these justified 
approximations, equation ( 2 ) becomes 

= (4) 

where /i e is the electronic contribution to the mean 
molecular weight. Then, dividing equation (4) by (nu), 
we arrive at the new equation 
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Moreover, taking the Taylor expansion around the av¬ 
erage temperature in phase i, (T) j, of the function of 
temperature f(T) given by the first term in claudators 
on the right-hand side of equation (5), we find that the 
average over ionized regions, i=II + III, of f(T) is well- 
approximated by /((T)i) + (d 2 //dT 2 )^ T ).cr^ i /2, where 
tTTi is the dispersion in temperatures around the mean. 

Besides being better justified than expression (3), ex¬ 
pression (5) is also more accurate for the following rea¬ 
sons: i) instead of taking the recombination rate den¬ 
sity at a fixed typical temperature divided by (/it®) hii, 
it uses the ^-dependent average of am (T)//i e in the 
Hii region, and ii) the last term on the right accounts 
for the changing comoving hydrogen density due to in¬ 
flows/outflows. In Figure 1, we compare (aHi(T)//x e )Hii 
for the ^-dependent temperature shown in Figure 2 to the 
uniform constant value aHi(Tty P )/(At e )Hii with T typ = 
10 4 K appearing in equation (3). As can be seen, the 
difference is noticeable, particularly around the redsliifts 
z = 10.3 and 5.5 of complete ionization in the particular 
galaxy model with double reionization considered. 

When Qhii reaches the value of one and (Ahii) is suf¬ 
ficient to balance recombinations, a period of ionization 
equilibrium begins in which Qhii stays equal to one. 



Figure 1. Average of the recombination to Hi coefficient over 
the electron contribution to the mean molecular weight, [i e , over 
Hii regions as a function of 2 for the temperature evolution shown 
in Figure 2 (solid red line), obtained from the galaxy model by 
Manrique et al. (2015) for realistic values of the parameters leading 
to double hydrogen reionization Salvador-Sole & Manrique (2014), 
compared to the usual value for a fixed temperature of 10 4 K 
(dashed blue line). 

However, if (Ahii) becomes insufficient to keep ionized 
regions growing (or stable), a recombination will begin. 
The constant and decreasing values of Qhii in those two 
regimes are also governed by equation (5), in the former 
case with (Ami) replaced by the equilibrium value, with 
the leftover metagalactic emissivity eventually used, duly 
redshifted, to ionize more hydrogen atoms at lower z’s. 

A similar derivation leads to the homologous equation 
for the Hem volume filling factor, 


dQHelll _ (AHelll) 
dt (n He ) 

aHeii(T) ^ 
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dt 
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Again, if at any point a period of Hem ionization equi¬ 
librium or recombination takes place, then Qiieiii stays 
equal to one or begins to diminish, respectively, accord¬ 
ing to the same equation ( 6 ). 


3. TEMPERATURE EQUATIONS 

Photo-ionization leads to photo-heating of the differ¬ 
ent IGM phases. Other heating mechanisms acting on 
the IGM are Compton heating by X-rays and by cosmic 
microwave background (CMB) photons at very high- 2 : 
(after decoupling of baryons from radiation at z ~ 150). 
Such heating is partially balanced by the cooling due 
to recombinations and desexcitations, cosmic expansion, 
Comptonization from CMB photons at low z, and col- 
lisional cooling (significant only in very hot neutral re¬ 
gions, if any). In addition, density fluctuations suffer 
gravitational contraction/expansion causing extra heat¬ 
ing/cooling. These are the main mechanisms causing the 
thermal evolution of the IGM. Below we mention (in ital¬ 
ics) a few additional mechanisms that are included in the 
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present more accurate treatment (see also Hui & Gnedin 
1997 for other possible heating and ionizing mechanisms, 
due to decaying or annihilating dark matter, not included 
herein). 

The local temperature of the IGM evolves according 
to the differential equation (e.g. Choudhury & Ferrara 
2005) 


dT 


dln(l + z) 
d ln/t 


= T 


2 dln(n/(n)j) 
3 dln(l + z) 
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dine 


din n 


d ln(l + z) dln(l + z) dln(l + Y)J 


( 7 ) 


The hrst term in claudators on the right, equal to 2, 
gives the cosmological adiabatic cooling of the gas ele¬ 
ment; the second term gives its adiabatic heating/cooling 
by gravitational compression/expansion for the baryon 
density n around the mean value (n); in region i, tak¬ 
ing into account that most diffuse IGM is in a linear or 
moderately non-linear regime; the third term gives the 
cooling due to the increase in mean molecular weight, 
p, caused by ionization and outflows from halos ; the 
fourth term gives the Compton cooling from CMB pho¬ 
tons, and the gain/loss of energy density, e, due to 
photo-ionization/recombination, Compton heating from 
X-rays, the achievement of energy equipartition by newly 
ionized/recombined fraction of gas (the different phases 
have distinct temperatures in general) plus mechanical 
heating accompanying outflows from halos ; and the fifth 
term gives the cooling/heating by the gain/loss of baryon 
density, n, due to outflows/inflows (this changes the av¬ 
erage specific energy of the IGM). As outflows take place 
from halos harboring luminous sources, we assume that 
they only affect ionized regions. 

Multiplying equation (7) by pn, and taking the aver¬ 
age over each specific phase under the approximation, for 
the reasons mentioned in Section 2, that p, e, n, and T 
do not correlate with each other, we arrive at 

dln ( T )i _ 2 dln((//)i(e)i/(n)i) 

dln(l + z) d ln(l + z) ’ 1 J 

with i=I, II, or III. Note that, in neutral regions (i=I), 
there are no stochastic effects of luminous sources: e does 
not change either through photo-ionization or by X-rays, 
p is kept strictly equal to the primordial value, and there 
is only a small change in n due to inflows. Consequently, 
a strong correlation is foreseen between the quantities e 
and n and temperature. Yet, we still ignore such a cor¬ 
relation for simplicity. This approximation is only nec¬ 
essary during the initial period of increasing ionization; 
in recombination periods, the gas properties in the new 
neutral phase remain uncorrelated as they have suffered 
important stochastic feedback effects from luminous ob¬ 
jects over the previous ionized phase. 

And what about the temperature dispersion around 
the mean in the different IGM phases, also required in 
equations (5) and ( 6 )? To calculate = ( T 2 )i — (T ) 2 
we need to consider the relation 


I dT2 _ T 2 L . 2 dln(n/(n)j) 

dln(l + ;r) _ 3 dln(l + z) 


d In (pe/n) 
dln(l + z) 


( 9 ) 


following from equation (7). The same steps above lead 



Figure 2. Average IGM temperatures in neutral (blue dotted 
lines), singly ionized (green solid lines), and doubly ionized (red 
dashed lines) regions obtained from the same model with double 
hydrogen reionization (at z = 10.3 and 5.5) and single helium 
reionization (at z = 2) as in Figure 1. Solid circles with error 
bars are the actual IGM temperatures estimated by Lidz et al. 
(2010) and (Bolton et al. 2010, 2012). 

to 


din (T 2 ) 1 / 2 _ 0 dln((/r)i(e)i/(n)i) _ dln(T)j 

dln(l + z) ^ dln(l + z) dln(l + *)‘ 1 ’ 

The initial conditions for equation (8) are (T)i (£;„;) = 
7cmb( -2-ini ) and (T)n( -2-ini) — (p)n{T)i( 

-2-ini )/(p) i, where 

z- ln [ is the redshift at which the IGM temperature be¬ 
gins to deviate from the temperature of CMB pho¬ 
tons 2 , satisfying 1 + 2 ini = 100(flb/i 2 /0.0125) 2 / 5 , 

where fib and h are the baryonic density parame¬ 
ter and the Hubble parameter scaled to 100 km s _1 
Mpc^ 1 . Similarly, the initial condition for equa¬ 
tion (10) is (T 2 )ii(z ini ) = (fj)ii{T 2 )i(zini)/{p)j , where 
(T^ 2 ) I (-2-ini) is equal to °tcMb(77j, -2-ini) *F "7q MB (2ini), 
where c4 CMB (z in i) = 27c MB (zini)<To(7?j, 2ii„i) 2 /3 is the 
CMB temperature variance at the Jeans scale at recom¬ 
bination, Rj, evolved to Zj n i , with ao(R,z) the 0-order 
spectral moment at the scale R and redshift z. The rea¬ 
son for the filtering at the scale Rj is that, at smaller 
scales, there were no temperature fluctuations at recom¬ 
bination, and the uniform temperature on those scales 
only suffered cosmological adiabatic cooling and the ef¬ 
fects of luminous sources, uncorrelated with T. If there 
is a period of increasing recombination, the initial mean 
temperature and variance in the recombined region are 
equal (except for different mean molecular weights) to 
those in the ionized phase giving it rise. We have checked 
that <r^((T)i) is always much less than (T) 2 , meaning 
that the second order Taylor expansion around (T); of 
any arbitrary function f of temperature is really close to 
the value f(T). 

In Figure 2, we show the temperature evolution 

2 Until that time, the residual density of free electrons and ions 
causes the gas to be thermalized by CMB photons. 
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that results from the present improved treatment of 
IGM evolution for the source functions, IV H i, iVnein, 
dln(e};/dln(l + z) and dln(n}idlii(l + z), provided by 
the same galaxy model as in Figure 1. 

4. HALO IONIZATION-BIAS 

The chance that halos with a given mass M at t 
will trap gas, and that the trapped gas will cool either 
through molecular bands or atomic lines and form metal- 
poor or metal-rich stars, respectively, depends on the 
temperature and ionization state of the IGM in which 
the halos are embedded. Consequently, the halo MF it¬ 
self must vary between neutral and ionized environments. 
Note that, given the homogeneity of the Universe, these 
probabilities are not a function of a specific point. In 
particular, the probability that a given arbitrary point 
lies in a ionized or neutral region is uniform and equal to 
Qffli(t) and 1 — QhiiW, respectively. 

To calculate the probability that a halo with mass M 
is located in an ionized region at the cosmic time t, 
P M (Hn,t), we will first consider the conditional prob¬ 
abilities Pm (H ii, i|H ii, tf) and Pm (Hii, f|Hi, tf) that the 
halo is in an ionized region at t given that it was either in 
an ionized or neutral region, respectively, at its formation 
at tf. The former of these two quantities is simply 

Pm( Hu, t|Hn, tf) = 1 - Pm(t, tf), (11) 

where Pui(t,t{) is the probability that the halo environ¬ 
ment recombines between tf and t because of the ab¬ 
sence of nearby sufficiently luminous sources. The latter 
is given by 

Pm{ Hii, t|Hi, tf) = Pm(^ tf) + Pau(t, tf), (12) 

where P^ (t, tf) is the probability that star formation be¬ 
gins to take place in a halo with M lying in a neutral 
environment between tf and t (we say “begins” because 
newborn stars soon ionize the medium around the halo), 
and Pan{t,t{) is the probability that the halo environ¬ 
ment will become ionized in the same period of time be¬ 
cause of the presence of nearby external ionizing sources. 

To derive equations (11) and (12) we have assumed 
that the probabilities Pm(Mf) and Pmi(Mf) are inde¬ 
pendent of halo mass. This may not be the case if there is 
some correlation between the halo mass- and ionization- 
biases. However, in terms of the effect of density on the 
ionization state of a region the tendency for halos har¬ 
boring more powerful ionizing sources to lie in higher- 
density regions contrasts with that for ionized bubbles 
to stretch more rapidly in lower-density regions, so they 
tend to balance one another. Therefore, even though 
the importance of this correlation is hard to assess with¬ 
out performing accurate hydrodynamic simulations with 
ionizing radiative transfer, we do not expect it to be too 
marked. In other words, the present treatment should 
be reasonably approximate. 

The total probability of finding a halo ionized at t can 
be expressed in terms of the above conditional probabil¬ 
ities and P M (Hn,t) upon formation, 

Pm (Hu, t) = P m (Hii, f|Hn, t f )P M (Hu, t f ) 

+P M (Hii,t|Hi,t f )[l-P M (Hii,tf)] . (13) 

Substituting the conditional probabilities on the right 
of equation (13) by expressions (11) and (12), setting 



Figure 3. Probability that star formation will begin to take place 
in halos with M at 2 : = 15 (blue lines on the right) and 30 (red 
lines on the left), before 100 Myr (solid curves), 30 Myr (dashed 
curves), and 10 Myr (dotted curves) after their formation in neutral 
regions. Note that at 2 : = 30 there are not yet any halos 100 Myr 
old. 



M (M e ) 

Figure 4. Halo MF in ionized environments (dashed line) and 
neutral ones (dotted line) at z = 15 (blue lines on the right) and 
30 (red lines on the left), compared to the global MF (solid line) at 
the same redshifts, obtained for the same realistic galaxy model as 
in previous Figures. The higher abundance in ionized environments 
of halos in a narrow range of low masses is due to the formation of 
new Population III stars in neutral regions. The rest of the halos 
in ionized regions arise from the ionization the previous objects 
produce around them. 


t = tf + At, and taking the limit of small At, equation 
(13) leads to the following differential equation governing 
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the evolution of Pm (Hn,i) 

dPM(Hn,t) dQmi(t) , t{) 

- Tt -= dt + d t [1 “ Qhii(<)] ' 

(14) 

To derive equation (14) we have taken Pm(t, tf) = 0 and 
Pmi(t, tf) = [Qmi(i) — <2ffli(if)]/[l — Qhii(^)] in peri¬ 
ods of increasing ionization, and Pm(t,t{) = [QHii(if) — 
and Tffli(Mf) = 0 , in periods of in¬ 
creasing recombination. Interestingly, in both cases one 
is led to the same differential equation (14), whose solu¬ 
tion for the initial condition Pm (Hii,0) = 0 yields the 
desired probability Pm (Hn,t) of finding a halo with M 
in a ionized region at t, its complementary value giving 
the probability Pm (Hi, t) of finding it in a neutral region. 

The probability -PJf(t,tf) in equation (14) is hard to 
estimate analytically because it depends on the number 
fraction of H 2 molecules, /h 2 , at the center of halos with 
M, whose PDF cannot be established without making 
appeal to the whole halo aggregation history. Thus, this 
function must be drawn from a full treatment of galaxy 
and IGM evolution. In Figure 3, we plot this function 
obtained from the same galaxy model as in previous Fig¬ 
ures. The halo MFs in ionized and neutral regions re¬ 
sulting from a global MF of the Slreth & Tormen (2002) 
form at two different redshifts are plotted in Figure 4. 
As can be seen, the higher the redshift, the more marked 
the effect,which is only visible, of course, before full ion¬ 
ization. 

5. SUMMARY 

In the present paper, we have derived an improved 
version of the master equations for the evolution of 
IGM ionization state and temperature, accounting for 
the composite, inhomogeneous, multiphase nature of this 
medium. Besides all the usual effects, the new version 
includes collisional cooling in hot neutral regions (nec¬ 
essary to deal with recombination periods as found in 
double reionization), mass exchanges between halos and 
IGM, and the achievement equipartition for newly ion¬ 
ized/recombined gas. In addition, we have derived the 
probability that a halo with a given mass M at 2 is 
located in a ionized or neutral environment, which is 
needed to accurately compute the source functions re¬ 
quired in the IGM master equations. 

To check the performance of this improved treatment 
of IGM we coupled it to the galaxy model by Manrique 
et al. (2015) for realistic values of the parameters leading 
to double reionization (Salvador-Sole & Manrique 2014). 
The main results were as follows: 

- The average temperatures in the three IGM phases 
show marked variations over the different ioniza¬ 
tion/recombination periods. This harbors relevant infor¬ 
mation on the epoch of reionization. The usual treatment 
dealing with the average temperature over the whole 
IGM (or at mean IGM density, To) loses this informa¬ 
tion. 

- The inclusion of collisional cooling is mandatory to re¬ 
cover the sudden decrement in the average temperature 
of neutral regions after first ionization in double reioniza¬ 
tion (see Fig. 2). In the only work to date, by Choudhury 
& Ferrara (2005), dealing with the evolution of the av¬ 
erage temperature in the different IGM phases, neutral 
regions cooled adiabatically after decoupling. 


- The average temperatures of singly and doubly ion¬ 
ized regions show a maximum similar to that found by 
Choudhury & Ferrara (2005; see panel f of their Fig. 1). 
However, our temperatures also show a minimum, due 
to the recombination after first ionization. More impor¬ 
tantly, the average temperature in doubly ionized regions 
is always higher than in singly ionized ones, while this 
was surprisingly not the case in Choudhury & Ferrara’s 
solution. 

- Although the average temperature in singly ionized re¬ 
gions is not as high as that reported by Choudhury and 
Ferrara, it is still notably higher (by a factor of ~ 3) than 
the value of 10 4 K often adopted in reionization studies 
(Shapiro et al. 1994; Wyithe & Loeb 2003; Benson et al. 
2006; Zhang et al. 2007). 

- This difference translates into the average recombina¬ 
tion coefficients. The values we find are substantially 
smaller (by a factor ~ 4) than found for the temperature 
of 10 10 K, and somewhat greater (by a factor ~ 2) than 
the minimum value at the average temperature reached 
in Choudhury & Ferrara’s solution. 

- This affects the evolution of the volume filling factors 
of ionized hydrogen and helium for identical source func¬ 
tions (identical galaxy models). But this makes a small 
difference compared to that arising from the galaxy mod¬ 
els used, which may lead, for instance, to single or double 
reionization. 

- We have computed the halo ionization-bias in the cal¬ 
culation of the source functions appearing in the IGM 
master equations. The ratios between the halo MF in 
ionized and all environments found for low mass newly 
star-forming halos and for the rest are respectively equal 
to ~ 3 x 1(T 4 (0.1) and ~ 0.3 (0.4) at 2 = 30 (z = 15). 

This improved treatment of IGM can be easily imple¬ 
mented in any model of galaxy and IGM evolution. This 
is particularly advisable for accurate models of galaxy 
formation or reionization when contrasting them with 
current observations (e.g. Salvador-Sole & Manrique 
2014) or future ones (e.g. 21 cm line experiments). 

This work was supported by the Spanish DGES grant 
AYA2012-39168-C03-02, and the Catalan DIUE grant 
2009SGR00217. 
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